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ABSTRACT 

By the virtues of the Dyson- Schwinger equations, we upgrade the published 
code HELAC to be capable to calculate the heavy quarkonium helicity ampli- 
tudes in the framework of NRQCD factorization, which we dub HELAC-Onia. 
We rewrote the original HELAC to make the new program be able to calculate 
helicity amplitudes of multi P-wave quarkonium states production at hadron 
colliders and electron-positron colliders by including new P-wave off-shell cur- 
rents. Therefore, besides the high efficiencies in computation of multi-leg pro- 
cesses within the Standard Model, HELAC-Onia is also sufficiently numerical 
stable in dealing with P-wave quarkonia (e.g. h Cjb , x c ,b ) and P-wave color-octet 
intermediate states. To the best of our knowledge, it is a first general-purpose 
automatic quarkonium matrix elements generator based on recursion relations 
on the market. 



PROGRAM SUMMARY 



Program title: 
HELAC-Onia. 

Catalogue number: 

Program obtainable from: 

|http: / /helac-phegas.web.cern.ch/helac-phegas 

Licensing provisions: none 

Operating system under which the program has been tested: 
Windows, Unix. 

Programming language: 
FORTRAN 90 

Keywords: 

quarkonium helicity amplitudes, NRQCD, Dyson- Schwinger equations, off-shell currents. 
Nature of physical problem: 

An important way to explore the law of the nature is to investigate the heavy quarkonium 
physics at B factories and hadron colliders. However, its production mechanism is still 
unclear, though NRQCD can explain its decay mechanism in a sufficiently satisfactory 
manner. The substantial K-factor in heavy quarkonium production processes also implies 
that the associated production of quarkonium and a relatively large number of particles 
may paly a crucial role in unveiling its production mechanism. 

Method of solution: 

A labor-saved and efficient way is to make the tedious amplitudes calculation automa- 
tion. Based on a recursive algorithm derived from the Dyson-Schwinger equations, the 
goal of automatic calculation of heavy quarkonium helicity amplitudes in NRQCD has 
been acheived. Inheriting from the virtues of the recursion relations with the lower compu- 
tational cost compared to the traditional Feynman-diagram based method, the multi-leg 
processes ( with or without multi-quarkonia up to P-wave states) at colliders are also 
accessible. 
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1 Introduction 



Studies of heavy-quarkonium systems,e.g. J /if), T and B c , provides an important opportu- 
nity to investigate quantum chromodynamics (QCD) at hadronic level with least artificial 
non-perturbative input parameters by hands. The fact relies on the non-relativistic prop- 
erty formed by relatively heavy charm and bottom quarks. Theoretically, these meson 
can be described well by non-relativistic QCD (NRQCD) [I] effective theory with only 
price that some non-perturbative parameters should be determined in prior. Although, 
in principle, the number of these parameters are not finite , for majority concerned phe- 
nomenology analysis, the number of important parameters are always limited given in 
velocity scaling rule. They can be determined once for all due to their universality prop- 
erty in the effective theory. In fact, NRQCD was established on a factorization theorem 
that a perturbative high-energy exchange part and a process-independent low-energy part 
works well in the production and decay processes of heavy quarkonium. The factoriza- 
tion has been proven rigorously for quarkonium decay to all orders, while the theorem in 
production processes is still absence of proof beyond two-loops. On the phenomenology 
side, the inconsistence of NRQCD factorization, cross section and polarization of heavy 
quarkonium production has not been eliminated yet [2J. Hence, fair to say, the mechanism 
of heavy quarkonium production is still unclear by now. 

The motivation of the paper is to develop an automatic tool for performing investiga- 
tions in heavy quarkonium physics. Compared to some of the traditional Feynman- diagram 
based tools[3J H],the package inherits the abilities of HELAC[5] El E], which is based on a 
recursive algorithm and hence reduces the computational cost that grows asymptotically 
as n! to a n with a ~ 3, where n is the number of external legs in the considered process. 
Therefore, we expect it provides a more efficient way for people to do physical anal- 
ysis with multi-particle processes, which might be especially important in quarkonium 
physics [8]. However, instead of making an extensive comparison between our package and 
others in the literature, we just want to focus on the realization of the automation with 
recursive relations in the framework of NRQCD. In section 2, the recursive algorithm 
in HELAC is revised. In section 3, we demonstrate the strategies in helicity amplitudes 
calculations of heavy quarkonium production in HELAC-Onia, and especially the descrip- 
tion of P-wave off-shell currents. Several benchmark processes are computed in section 4. 
Finally, we explain our program and draw our conclusions and outlooks in the last two 
sections respectively. 
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2 The recursive algorithm 



The algorithm of HELACfSJEJE] is based on the Dyson- Schwinger equations JHJQIJJQT], which 
is an generalized version of the Berends-Giele off-shell recursive relation [T2]. To illustrate 
it, we consider a process with n external legs. The momentums of these external legs are 
denoted as pi,P2, ■ ■ ■ ,Pn, and the quantum numbers (e.g. color, helicity) are defined as 
cki, «2> • • • j o: n . An off-shell current with k external legs can be represented as 



J({p il ,...,p ik };{a il ,...,a ik }) 




P = Pii + Pi-2 + • • • + Pi k 



Pi k , a i k 



All of the subgraphs that are able to transfer the k external legs into the off-shell current 
J according to the Feynman rules in the considered model have been included into the 
shade bubble. Every current is assigned by a "level" I, which is defined as the number 
of external legs involved in the current, i.e. the "level" of ^({p^, ■ ■ ■ ,Pi k }', {a^, . . . , a ik \) 
is k. In special, the "level" of each external leg is 1. Then, in general, all of the currents 
with higher "level" can be constructed from those with lower "level". The starting point 
of the recursion relation is from the / = 1 currents, i.e. the external legs. For the i-th 
external leg, the corresponding currentQ is its wavefunction 

J({PiY, {at}) = > (2) 

Specificallyfor a vector boson 

JiWiinAY) = (3) 

where \i is the lorentz index and A is the helicity of the vector boson (A = ±1 for a 
massless vector, while A = ±1,0 for a massive vector )0, while for a spin-| fermion 



J({ Pi };{+l,X}) 



u \(.Pi) when p° > 
V\(—Pi) when p° < 



J(M;{-i,A}H !\ W , , (4) 

I V\{— pi) when pj < 



1 The "level" 1 current is on-shell instead of off-shell. 

2 Note that, for simplicity, we have suppressed other possible quantum numbers like color for gluon. 
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where +1 and — 1 means fermion flow and anti-fermion flow respectively, A is the helicity 
index with A = ±1. The explicit expressions of these / = 1 wavefunctions are presented 
in the appendix of Ref.|5]. The currents with I — k > 1 can be constructed the currents 
with lower |] 




where a means to exhaust all possible generating " level" r, s (and t) currents formed by 
the ii, . . . ,ik external legs. Each off-shell currents should be multiplied by its propagator. 
The end of the recursion is the forming of the "level" n current, in which we choose all 
I — n— 1 currents^ to multiply with the first external particle's wavefunction. If the flavor 
of the first external particle is not exactly same as the flavor of the I = n — 1 current, 
the current is dropped. Finally, we obtain the resulting amplitude. In this way, one 
can avoid computing identical sub-graphs contributing different Feynman diagrams more 
than once. The summation of the all sub-graphs that contribute to a specific current also 
reduces the number of objects that should be used in the next "level" recursion formula. 
Therefore, the computation complexity will be reduced from ~ n! in Feynman- diagram 
based algorithm to ~ a n in the Dyson-Schwinger based recursive algorithm, where a ~ 3. 

In the original HELAC program [5], it uses a binary representation of the momenta in- 
volved in the considered process[l3]. For each of the external momenta Pi, ■ ■ ■ ,p n , its 
binary representation is 2 l ~ l for z-th external leg with momenta Pi, while for a I = k 
current J7"({Ph, ■ ■ ■ ,Pi k }] {ot^, ■ ■ ■ , a ik }), it is expressed as J2j=i 2 lj_1 . Then, every mo- 
menta P M = J2]=i m j Pj can be uniquely expressed by an integer m = J2]=i 2 J-1 m ; / where 
rrij = or 1. In this case, the "level" of an current with momenta P M = Y%=i m j Pj can be 
calculated directly by I = J2]=i m j- 111 this case, the sign factor from the anti-symmetric 

3 For simplicity, we only consider tri-linear and quadri-linear couplings. However, it is straightforward 
to include higher-point vertices as well. 

4 Actually, they are on-shell instead of off-shell. 
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property of fermions is obtained by 



-i 



with 



e(P 1 ,P 2 ) = (-l)^ P2 ), X (Pi,P 2 ) = E^E^i (6) 

i=n j=l 



3=1 



P 2 = J2 m 2jPj> 
3=1 

{0 when particle j is a boson 
mij when particle j is a fermion 



my 



A | when particle j is a boson 

2j | when particle j is a fermion 

If the current is constructed by a tri-linear coupling with the lower " level" currents Pi and 
P 2 , it should be multiplied by a factor e(Pi,P 2 ). If it is constructed by a quadri-linear 
coupling with currents Pi,P 2 and P3, the sign factor is e(Pi,P 2 ,Ps) = e(Pi,P 2 )e(Pi + 

The way of the color treatment is also an interesting topic in the matrix element 
generator. In HELAC, it is using the widely used color flow basis, which was first proposed 
in Ref. [H] and was applied in perturbative QCD computations in Refs. [T5l [16] . Basically, 
a color octet gluon filed A® is replaced by a 3 x 3 matrix (-4 M )* = ^A^(A a )*-, where A a is 
the Gell-Mann matrix, i.e. 8 = 3 £g> 3 — 1, while incoming quarks or outgoing antiquarks 
still maintain in the 3 representation of SU(3) and outgoing quarks or incoming antiquarks 
are in 3 representation. After this substitution, only the Kronecker notation 5s appear in 
the Feynman rules. All the Feynman rules in the color-flow basis have been established 
in Ref. [16]. If there are n g external gluons (denote as 1, 2, . . . , n g ) and n q external quark- 
antiquark pairs (denote as n g + 1, n g + 2, . . . , n g + n q ) in the considered process, the color 
basis for the amplitude will be in the form of 

C i = ■ ■ ■ ^K+r^)' ( 8 ) 

where ex, represents the i-th permutation of 1, 2, . . . , n g + n q . There are totally {n g + n q )\ 
color basis ,though some of them will vanish. With this basis, one can construct the color 
matrix via 

m^ = 5>c*, (9) 
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and obtain the final square of matrix elements by 

(n g +n q )\ 

\M\ 2 = £ AiMijA*, (10) 

where A4 and Aj are the color-stripped amplitudes. 

In order to improve the computation efficiency, a Monte Carlo sampling over the he- 
licity configurations is adopted in the program [5] to perform the helicity summation. The 
basic idea of this technology is simple. Let's take a massive vector boson for example. 
A massive vector boson has three helicity states A = ±1,0 with wavefunctions e+,e^,eg. 
The strategy puts the concrete helicity summation of J2x=±,o e \( e \)* a continue in- 
tegration by defining = Sa=±,o e tX ^e^. Then, the summation becomes Jq 77 d^e^(e^)*, 
which can be calculated easily by a Monte Carlo program. 

3 Quarkonium amplitudes in NRQCD 

In the framework of the NRQCD factorization, the cross section of a heavy quarkonium 
production can be written as a combination of the perturbative short- distance parts and 
the non-perturbative long-distance matrix elements. For example at the proton-proton 
collider, the factorized form of a heavy quarkonium Q production is written as 

a(pp^Q + X) = Y,J dx 1 dx 2 / i/p (x 1 )^ /p (x 2 )a(ij ^QQ[n]+X)(0 n s >, (11) 

where f^ p and fj/ p are the parton distribution f unctions, a (ij — > QQ[n] + X) is the short 
distance cross section of producing a heavy quark pair QQ in a specific quantum state 
n, and (O^) represents as the long distance matrix element. In principle, for a specific 
quarkonium Q, there are infinity number of Fock states n and infinity number of long 
distance matrix elements (Of~). The power counting rules in NRQCD tell us for any 
quarkonium, there are only limited Fock states should be involved in our calculations up to 
a specific order of v , where v is the relative velocity of the heavy quark pair that forms the 
quarkonium. It makes NRQCD be predictable for hadrons. For example, in the process 
of J/ip production, there are only four different Fock states (i.e. B, ^Sf\ 3 Pj 8 ^ and ^q 8 '^ 
contributing to its cross section up to v 7 . The color-singlet long distance matrix element 

5 We write the Fock states in the spectroscopic form of n — 2S+1 L^j\ where S, L, J identify the spin, 
orbital momentum, total angular momentum states respectivley, and c = 1, 8 means that the intermediate 
state QQ can be in color-singlet or color-octet state. 
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can be estimated from the phenomenological models like potential models, while the color- 
octet long distance matrix elements can only be determined from the experimental data 
till now. 



3.1 Projection method 

To evaluate the pro cess- dependent short distance coefficients, one has to constraint the 
QQ into a specific quantum state. A convenient way to do it is performing projection. 

Specifically, the color projectors to the process ij — > QQ[ 2S+1 L^-j] +X are [IT] jf- when 
c = 1 and y/Ztfj when c = 8, where i, j are the color indices of the heavy quark pair QQ 
and A a is the Gell-Mann matrix. The color octet projector which contains the Gell-Mann 
matrix will be decomposed into the color-flow basis in the HELAC-Onia. Moreover, after 
projecting, no color indices of the heavy quark pair in the color-singlet states will appear. 

Another important constraint of the heavy quark pair is their total spin. The spin 
projectors were first derived in Ref.[18]. The general form of the projectors is^| 

-v(p 2 ,X 2 )T s ^ n^,^), (12) 



2V2{E + m Q ) y ' 2E 

where rriQ is the mass of the heavy quark,pi,p 2 and Aj, A2 are the heavy quarks' momenta 
and helicity respectively, P^ = Pi + p% is the total momentum of the heavy quark pair 
and E = The Y $ is 75 when S — 0, and it is e^ s 7 M when S = 1, where A s = ±, 

is the helicity of the quarkonium Q and e^ s is the polarization vector for the spin triplet 
state. For S-wave and P-wave states without relativistic corrections, E can be set as 
rriq directly. After applying the spin projection, the two external wavefunctions for open 
Q and Q will be glued. It results in a problem in the recursive relation, because the 
recursion begins from the external wavefunctions. In order to cure it, we decide to cut 
the glued fermion chain at the place of f + 2E in the projector Eq. (|12p . Using the 
completeness relation of f + 2E = J2y=±u(P, X')u(P, A'), we use the new "wavefunction" 
for Q as ^u(P, X')(f l + itlq) and for Q as — s ^ mQ {f 2 — mo)u{P, A'). Considering the A' 
in the "wavefunctions" of Q and Q should be exactly same, we have to perform the direct 
summation of A' in stead of Monte Carlo sampling in the HELAC-Onia. 



3.2 P-wave currents in HELAC-Onia 

The P-wave calculations are always necessary in the NRQCD predictions for both P- 
wave states h c /b,x c /b and S-wave states J/ip,T,i] c /b. For example, the color-octet P- 

6 In HELAC-Onia, we also generalize the projectors in the case of the heavy quarks in different flavors 
that form a heavy quarkonium like B c . But for simplicity, we only consider the same flavor case here. 
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wave states 3 Pj are playing special roles in J/ip hadroproduction [T9"} I2TI I23] and 
photoproduction[24, 25J. Although they are power suppressed in NRQCD compared to fi, 
fragmentation topologies make them overwhelming the color-singlet one at the medium 
and high transverse momentum regime. Hence, HELAC-Onia is designed to be able to 
handle with P-wave states as well with a numerical stable method by introducing new 
P-wave off- shell currents. 

After expanding the relative momentum q v = 1 2 between the heavy quark pair in 
the amplitude A(ij — > Q(pi)Q(p2)+X) in the non-relativistic approximation, the formula 
for the calculation of P-wave amplitude is 



where A/ = ±, is the helicity configuration of the polarization vector for P-wave state. 
The treatment of the new "wavefunctions" definition of the heavy quark pair avoids the 
derivation of the spinors, which might result in numerical instability. 

Alternatively, one could also do a direct numerical derivation by keeping the small 
relative momentum q of the quark and anti-quark that forms the heavy quarkonium and 
approaching q to zero in the quarkonium rest frame [3j. However, this direct numerical 
derivation might result in numerical unstable potentially especially when there are many 
P-wave states involved in the process. 

In contrast, the P-wave currents, which are extended from the original off-shell currents 
at parton level, can be written in a much compact manner. In the HELAC-Onia, we assign 
each current with an derivation index, which is also in binary representation. Assuming 
there are np P-wave states in the considered process, the relative momenta of the z-th 
heavy quark pair that forms P-wave state is denoted as q,i where i = 1, . . . , rip. The 
general derivation index form for a current is b = Y^i=i bii % ~ 1 with bi = or 1. If the 
current has been derived by qi as done like in Eq. (iT5|) .&,- is 1, otherwise bi = 0. Finally 
, only the amplitudes with b = 2 np — 1 are kept. The numerical stable form of P-wave 
currents avoids the large numerical cancellation. 

4 Benchmark processes 

We are in the position to illustrate the validation and applications of HELAC-Onia to 
the heavy quarkonium production at the proton-proton, proton-antiproton and electron- 
positron colliders. 




_d_ 

dq, 



V 



A(ij^Q(p 1 )Q(p 2 )+X) 



(13) 
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4.1 B c meson production at the LHC 

B c production is an interesting channel to investigate QCD, and it has been widely studied 
at the hadron colliders [26j ETJEE]. The available results have been used by the MAD0NIA[3] 
for testing the correctness of the code. We will also compare our results calculated by 
the HELAC-Onia with those presented in Ref. [3]. We only consider the B c production at 
the LHC with the center-of-mass energy 14 TeV with the initial gluon-gluon fusion and 
quark- ant iquark anihilation here. All of the input parameters are taken as same as those 
in Ref. [3]: 

(1) The masses of the bottom, charm quarks and B c meson are set as m& = 4.9 GeV, m c = 
1.5 GeV, m Bc = m b + m c = 6.4 GeV. 

(2) The parton distribution function (PDF) set is chosen as CTEQ6L1 [29J. 

(3) The factorization scale //p of PDF and the renormalization scale (Ir are set as 
/J-f = /J-r = 2(m{, + m c ) = 12.8 GeV. Moreover, the strong coupling constant is fixed 
as ois(hr) = 0.189. 

(4) The color-singlet long distance matrix elements are taken as (O b °( 2S+ *S [ } ] )) = (2J + 
1)0.736 GeV 3 , (0 B °( 2S+i P l j ] )) = (2J + 1)0.287 GeV 5 . 

(5) The color-octet long distance matrix elements are related with the color-singlet 
ones,i.e.(0^( 25+1 5' 8] )) = (O^+^/lOO, (O b *( 2S+1 P? ] )) = (0 B <f s+1 pj 1] ))/100. 

Our final results (with Monte Carlo statistical errors) are shown in the second column of 
TableJU where we also listed the corresponding results (in the third column of TableJT)) 
presented in Ref. [3] for the convenience of comparison. We find that our results are in 
agreement with those in Ref. [3J. 

4.2 Charmonia production at the B factory 

The charmonia production from the electron-positron collisions has been extensively stud- 
ied over the past decade. We don't intend to recall the long story of the theoretical and 
experimental studies on this topic here, which was already summarized in Ref.|30j. We 
only want to show the application and validation of the program HELAC-Onia in calculat- 
ing the quarkonium observables in the e + e~ colliders in this section. 

The first theoretical results of the inclusive charmonium association production with 
cc via single virtual photon exchanging at the B factory with the center-of-mass energy 
10.6 GeV was presented in Ref. |31j, which are only calculated at the leading order in as 
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process 


EhLAC-Oma(nbJ 


Tiff A Pv ATVTT A ( 1 \ 

MADONlA(nbJ 


99 


— y 




)oc 


39.3994 ± 0.0958382 


OA A 

39.4 


99 


— y 




\ I,— 

joe 


98.3109 ± 0.287252 


98.3 


99 


— >• 


m-4- /In [11 


\ 7 — 

Joe 


5.21131 ± 0.0144431 


5.20 


99 


— y 


B c( P J 


Joe 


16.7341 ± 0.0589108 


16.72 


99 


— y 




)oc 


0.411671 ± 0.00169734 


0.411 


99 


— >■ 




\ 7. — 

)oc 


1.78657 ± 0.00624756 


1.79 


99 


— y 




Joe 


0.11816 ± 0.000754526 


0.117 


99 


— y 


7-)_i_ /3t-)[8] 


\ 7 — 

Joe 


0.305862 ± 0.0011841 


0.3051 


qq 


— y 


B c{So 


\ 7 — 

)oc 


0.137782 ± 0.000896985 


0.137 


qq 


—y 


B c 1^1 


)bc 


0.83905 ± 0.00524885 


0.834 


qq 


— y 




)bc 


U-Uzyolzo ± U.UUUlo4yiy 


U.Uzyo 


qq 


->• 




)6c 


0.111259 ±0.000839535 


0.1105 


qq 


— >■ 


B tes [ o ] 


)6c 


0.00103294 ±4.44716 • 10~ 6 


0.00103 


qq 


— >■ 




)6c 


0.00707624 ± 0.0000459292 


0.00703 


qq 


— >■ 




)6c 


0.000253678 ± 2.19206 • 10~ 6 


0.000251 


qq 


->■ 


£ c+( 3p[8] 


)6c 


0.000826534 ± 5.16988 ■ 10" 6 


0.0008207 



Table 1: Cross sections of inclusive B+ production at the LHC with the center-of-mass 
energy 14 TeV. The data in the third column are taken from Ref.[3]. In the second column, 
the Monte Carlo statistical errors are also given. 

and v. Moreover,the results of rj c and J/ip production with gluons at the B factory are 
given in Ref.[3]. We put the same input parameters into HELAC-Onia: 

(1) The charm quark's mass m c is 1.5 GeV, and the masses of the charmonia considered 
are 2m c = 3 GeV. The mass of the electron and positron are safely ignored. 

(2) The renormalization scale fiR is chosen as 2m c = 3 GeV. In this way, the strong 
coupling constant is fixed as as((iR) = 0.26, while the electromagnetic fine structure 
constant is also set as a = 1/137. 

(3) The color-singlet long distance matrix elements are (C( 25+1 S'j')) = (2 J+l)0.387 GeV 3 . 
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Our color-singlet S-wave results are presented in Tablej2j from which we see that all of 
our results are in good agreement with those in Refs. [3TJ [3] . In first two rows of TableJU 
we also presented the cross sections of r\ c and J/ip production in association with cc at 
0(a 2 a'g+a 3 as+<y 4: ), which means that we have include both the single photon exchanging 
channel and the double photon exchanging channel. To the best of our knowledge, they 
are new. 

Several years ago, it was reported that there was a large discrepancy between the theo- 
retical predictions and the experimental measurements in the exclusive double charmonia 
production at the B factory (see review in e.g. Ref. [30]). The large discrepancy has at- 
tracted a lot of studies in this fields. We recalculated some of the cross sections at the B 
factory with yfs = 10.6 GeV in Tablej3]with the same input as those given in Ref. [32]. 
We listed the input parameters as following: 

(1) The mass of the charm quark is 1.4 GeV, and the masses of the charmonia are 2m c . 

(2) The renormalization scale hr is set as \fs/2 = 5.3 GeV, and asi^n) = 0.21, a = 
1/137. 

(3) The color-singlet long distance matrix elements are (0( 2S+1 S [ j ] )} = (2 J + 
1)0.335 GeV 3 and (C( 25+1 Pj 1] )) = (2 J + 1)0.053 GeV 5 . 

The results in Tablefjonly include 0(ot 2 a 2 s ) and leading-order in v perturbative calcula- 
tions. Good agreement is found between HELAC-Onia and Ref.|32j. The 0{a 2 a 2 s + a i ots + 
a 4 ) cross sections for J/ifjJ/ifj and J/iph c exclusive productions are shown in the last two 
rows of TableJU These results agree with those given in Ref. |33j . 



process 


HELAC-Onia(fb) 


Refs. [311 E] (fb) 


e+e" ->• 7* ->• VcCS l o)cc 
e + e~ -»■ 7* T] c CS [ o ] )ggg 
e+e" -> 7* -»■ J/i)(*S [ l ] )cc 
e+e" -> J/tPCS [ i ] )99 


58.7938 ±0.154193 
3.72893 ±0.0063512 
147.864 ± 0.305001 
266.037 ±0.247366 


58.7 
3.72 
148 
266 



Table 2: .Cross sections of the inclusive r] c and J /if} production via single virtual photon 
exchanging at the B factory with the center-of-mass energy 10.6 GeV. 
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process 


HELAC-Onia(fb) 


Ref. 


[32] (fb) 


e+e - _> Y ^ ( ip[i])^(i5ji]) 
e+e - ^ 7 * ^ JM^Xcjfif) 
e+e - _> 7 * _> ^pf^fpj 11 ) 


3.78154 ±0.00338108 
0.308533 ±0.000198459 

3.47635 ± 0.00453553 
0.328299 ±0.000392734 


3.78 
0.308 
3.47 
0.328 



Table 3: .Cross sections of the exclusive double charmonia production via single virtual 
photon exchanging at the B factory with the center-of-mass energy 10.6 GeV. 



process 


HELAC-Onia(fb) 


Ref.[33j(fb) 


e+e" ^ Vc Cs£ l )cc 
e+e- -> J/if)(*S [ p)cc 
e+e- -> J/ij( 3 S\ L] )J/ij( 3 S\ L] ) 
e+e- -> J/^(^)/ lc ( 1 P 1 W) 


61.6802 ±0.0854359 
166.499 ±0.175318 
6.64805 ± 0.0123474 
0.00606923 ± 6.84416 • 10" 6 


6.65 
0.0061 



Table 4: .The 0{a 2 a 2 s + a 3 as + a 4 ) cross sections of the charmonia production at the B 
factory with the center-of-mass energy 1.96 GeV. 

4.3 Double quarkonia production at the Tevatron and the LHC 

Double quarkonia production at the hadron colliders is a useful way to investigate the 
color-octet mechanism. In this section, we will compare the results calculated by HELAC-Onia 
and those in the literature [34]. The input parameters (same as those in Ref. [33]) are: 

(1) m c = 1.5,m6 = 4.9. The mass of a heavy quarkonium is just approximated as the 
sum of its constituent heavy quarks' masses. In other words, if a heavy quarkonium 
H is composed by Q1Q2, then uih = r mQ 1 + mQ 2 - 

(2) up = fM R = y m 2 H ± p\ for the heavy quarkonium H . 

(3) PDF set is CTEQ6L1 [29]. Therefore, the running of as is evaluated by the leading- 
order formula in the PDF set. 

(4) The S-wave color-singlet long distance matrix elements are 

1)0.389134 GeV^O^+^i 11 )) = (2J±1)2.34722 GeV 3 and (O b5 ( 2S+1 S]j ] )} = (2J+ 
1)0.720017 GeV 3 . 

(5) The pseudorapidity rj cuts on the final quarkonia are \rj\ < 0.6 at the Tevatron and 
\rj\ < 2.4 at the LHC. 
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Final States 



2T( 3 Sj 1] ) 

B e frjp)B c (*5P) 
5 C (M 1] )5 C (M 1] ) 



HELAC-Onia(nb) 



Ref.[34j(nb) 



3.32 ■ IO" 3 

0.0563 
1.87- 10~ 5 



3.316- 10~ 3 ± 3.705- 10" 
0.05631 ±4.437- 10~ 5 



1.866 
1.226 
3.854 
1.001 
8.226 



IO" 5 ±2.385 
IO" 4 ± 1.489 
10" 3 ± 9.529 
10~ 3 ± 2.492 
10~ 3 ± 9.531 



io- 8 
io- 7 

IO" 6 
10~ 6 
10~ 6 



1.23 
3.86 
1.00 
8.23 



io- 4 

IO" 3 
IO" 3 
10~ 3 



Table 5: Cross sections of double quarkonium production at the Tevatron with the center- 
of-mass energy 1.96 TeV. 



Final States 


HELAC-Onia(nb) 


Ref.l34j(nb) 


2JM^1 1] ) 

2%( 1 4 11 ) 

2Y( 3 S , {^) 
B^BcCsP) 


2.730 ±0.01710 
2.832 ± 1.721 ■ IO" 3 
7.373 • 10~ 3 ± 1.802 • 10~ 5 
0.01514 ± 1.184 • 10~ 5 
0.2723 ± 1.461 • IO" 4 
0.08379 ±4.430 • 10~ 5 
0.7078 ± 3.797 ■ IO" 4 


2.73 

2.83 
7.36 • 10~ 3 
0.0151 

0.272 
0.0837 

0.708 



Table 6: Cross sections of double quarkonium production at the LHC with the center-of- 
mass energy 14 TeV. 

The S-wave color-singlet cross sections are shown in Table J5] (Tevatron with y/s = 1.96 TeV) 
and in TableE] (LHC with y/s = 14 TeV). 



4.4 Hadroproduction of J/ip and T in association with a heavy- 
quark pair 

The measurements of the J j ip and T in association with a heavy quark pair at the hadron 
collider are interesting because not only they will contribute to the inclusive J ftp and T 
production but also they are useful way to study color-octet mechanism at the Tevatron 
and the LHC. In figure [U, we present the transverse momentum distributions of J/ip 
and T via color-singlet channel at the Tevatron with \fs = 1.96 TeV and the LHC with 
y/s = 14 TeV. All of our results are in agreement with those in Ref . [35] . For completeness, 
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we list our input for calculation J/ipcc and Tbb production as follows: 
(1) m c = 1.5 GeV,nib = 4.75 GeV and mj^ = 2m c , my = 2m&. 



(2) = = y (4wiq) 2 + p?p } where uiq is m c for J/ipcc and m& for Tbb. 

(3) PDF set is CTEQ6M [29]. The running of as is following the next-to-leading order 
formula in CTEQ. 

(4) The color-singlet long distance matrix elements are (O j ^(Q)) = 1.16024 GeV 3 and 
(£> T (B)) = 9.28192 GeV 3 . 

(5) The rapidity cuts are applied as \y\ < 0.6 at the Tevatron and \y\ < 0.5 at the LHC. 



4.5 Spin density matrix and polarization 

Besides the total cross sections and other unpolarized observables like pt spectrum, 
HELAC-Onia is also designed to be able to calculate the spin density matrices of heavy 
quarkonia. Hence, it can be taken as a useful tool to calculate the polarization observables 
of heavy quarkonia in various polarization frames. It has been successfully used in: 

(1) the next-to-leading order inclusive J/ip polarization at the Tevatron and LHC [22J; 

(2) the polarization of the inclusive \c hadropro duct ion [36| 157]; 

(3) the polarized \ c production in associate with a charm quark pair at the LHC|38]. 

The readers who are interested in these topics can refer to the corresponding (forthcoming) 
publications. 

Besides the above examples, there are many other aspects of quarkonium physics one 
can perform analysis with HELAC-Onia, for instance, the quarkonium production with jets 
(inclusive or exclusive^) and/or weak bosons. We refrain to illustrate more examples here. 

One can feel free to use the modified PHEGAS [HI [7J,RAMB0 [42] or VEGAS [43] to perform 
Monte Carlo evaluation. Standard Les Houches Event files [44] are also generated. 

Thanks to the useful discussion with Qiang Li, we have implemented the matrix element and parton 
shower matching method MLM scheme [39, 40 in the program. 



15 



10 2 
10 

1. 

2 10-1 

1 io- 2 

b 

X! 1Q -3 

io- 4 

IO -5 



J I if/ c c 
Color Singlet 

LHC 

Tevatron 








10" 



t 

Q. 

-a 

b 



10" 



10" 







1 1 






r ~- 

— \ 






Xhb 








Color Singlet 






\ 

\ 




LHC 






\ 


"S. 










\ 

\ 


Tevatron 










\ 

N, 










\ 

N 

\ 

\ 










\ 

\ 










\ 

s 















10 



15 

Pt (GeV) 



20 



25 



30 



Figure 1: pt distributions of J/tp(T) production in association with cc(bb) at the Tevatron 
and LHC. Only color-singlet states are considered. 
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5 Running the program 



We are in the position to describe how can one run the program. The program is 
split in two major phases, which were called initialization phase and computation phase 
as in Ref. [5]. During the initialization phase, the program selects all the relevant sub- 
amplitudes for the required process and evaluates the color matrix My, while during the 
the computation phase,it computes the amplitude for each phase space point introduced 
by PHEGAS,RAMBO or VEGAS. 

The running of the program is very easy. If the program is running under the Unix, 
one should specify the Fortran90 compilation in the first line of makefile. The default 
one is gf ortran. For the Windows user, it is running only after the user has included all 
of the Fortran90 files in his/her project. There are two input files that should be specified 
by the user before running the program. They are process . inp and user . inp. 

In the process . inp, the user should tell the program the information of the process 
including the number of external particles (in the first line) and the ids of the particles 
(in the second line). The table of the ids for the particles ( not hadrons ) in the standard- 
model defined in HELAC-Onia are the same as that in HELAC, which is shown in Tabled 
The naming rules of the ids for the heavy quarkonia in HELAC-Onia are: 

(1) The ids of the heavy quarkonia are in 6-digits. 

(2) The first two digits are 44 for charmonia, 55 for bottomonia and 45 for B c . 

(3) The next four digits are just recording the information of which intermediate Fock 
states. In general, the four digits are in the order of (2S + l)LJc for 2S+1 L [ } . For 
example, 3118 means the intermediate state is P[ . 

(4) Charmonia and bottomonia are all self-conjugated mesons, while B c are not. A 
minus sign is used to represent the anti-particle. In the program, we treat 

o ro] 

as the particle while B~ as the anti-particle. For example, the id of B~(P{ J ) is 
-453118. 

Using these rules, one can calculate the helicity amplitudes for S-wave and P-wave quarko- 
nia production from pp,pp and e + e~ collisions. We take one example. If the user want to 

q fa] 

calculate gg — > cc[P x ] + cc. The first line of process, inp is 5, and the second line of 
process . inp is 

35 35 443118 7-7. 
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u e ,e ,u,d,u^,fi , c, s, ... 


1,.. 


.,12 


i7 e ,e + ,u,d,Vn,fi + ,c, s, ... 


-I,-- 


.,-12 


"f,Z,W+,W-,g 


31,. 


..,35 


H, X, + , 0" 


41,. 


..,44 



Table 7: The ids of the "elementary" particles in the standard model in HELAC-Onia. 

The file user . inp is left for the user to specify the parameters in default . inp if he/she 
doesn't want to use the default values given in default . inp0. The main parameters are: 

(1) colpar represents for the type of colliding particles, i.e. 1 for pp,2 for pp and 3 for 

(2) energy is the center-of-mass energy y/s in unit of GeV. 

(3) gener specifies the Monte Carlo generator, i.e. for PHEGAS, 1 for RAMBO, 2 for 
DURHAM and 3 for VEGAS. 

(4) ranhel is a parameter to determine whether the program uses the Monte Carlo sam- 
pling over the helicity configurations. In specific, if ranhel= 0, it does the helicity 
summation, while if ranhel> 0, it does the Monte Carlo sampling. If ranhel= 1, 
the program uses Monte Carlo sampling over the helicities of the "elementary" 
particles in the standard-model and summing over helicities of quarkonia, while if 
ranhel= 2 it also performs Monte Carlo sampling over e£ for the P-wave states, 
and ranhel= 3 means it does Monte Carlo sampling over all polarization vectors of 
heavy quarkonia (of course also over helicites of the "elementary" particles in the 
standard- model) . 

(5) The value of qcd determines the amplitudes should be calculated in which theoryi.e. 
for only electroweak, 1 for electroweak and QCD, 2 for only QCD, 3 for only QED 
and 4 for QCD and QED. 

(6) alphasrun is a parameter to determine whether the strong coupling constant as 
should be running(l) or not(0). 

(7) Flags like gauge, ihiggs and widsch determine the gauge(0 = Feynman gauge, 1 = 
unitary gauge), whether inclusion Higgs(l) or not(0) and using the fixed(O) or com- 
plex^) scheme for the widths of W ± and Z bosons. 



8 We suggest the user don't change the content in the file default, inp unless he/she really knows 
what he/she is doing. 
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(8) nmc is the number of the Monte Carlo iterations. 

(9) pdf is the PDF set number proposed in LHAPDF[35]. Entering means no PDF 
is convoluted. 

(10) ptdisQ is a flag whether the pt distribution of the first final quarkonium are calcu- 
lated(T) or just total cross section(F). If ptdisQ is T, one should also specify which 
Pt value (Ptl) should be calculated. 

(11) Scale specifies which renormalization (and PDF factorization) scale should be used. 
It is explained in the comment line of default . inp. If the user chooses the fixed- 
value scheme, he/she should also supply the value of the scale (FScaleValue). 

(12) exp3pjQ is a flag whether summing over(F) Pj , J — 0, 1, 2 or not(T). 

(13) modes determines whether the calculated result is the polarized one(l) or not(0). If 
it is 1, the user should also supply the values of SDME1 and SDME2 to let the program 
know which spin density matrix element to calculate. Meanwhile, the value of LSJ 
represents which "spin" in quarkonium should be specified. The user should also 
specify the polarization frame (PolarFrame). 

(14) The parameters of the physical cuts in calculating the cross sections are also should 
be input by the user if he/she wishes to use his/her values. 

(15) The long distance matrix elements are also supplied in default, inp. The user 
can supply his/her values in user, inp with the same format in default . inp. The 
conventions are explained in the comment lines of default . inp. 

All other parameters are listed in default, inp. The user can fix his/her values in 
user, inp following the format in default, inp. Finally, the user just run the program 
and obtain the result. 

6 Summary and outlooks 

The exploitation of the fundamental law in the nature is the main aim of the running of 
the LHC Heavy quarkonium, which is one type of the simplest hadrons, provides an ideal 
laboratory to test and understand QCD. The discrepancies between the experimental 
measurements and the theoretical predictions imply that we still don't understand the 
production mechanism of the heavy quarkonium. Moreover, the quarkonia like J / ip and T 
production at the LHC will not only be taken as calibration tools but also be very useful 
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for TeV physics and even new physics. Hence, it is mandatory to improve the reliability 
of Monte Carlo simulation. From the present theoretical studies, the radiative corrections 
are very indispensable even at qualitative level to the quarkonium production. 

Unlike the case in parton-level, the automatic tools for calculating quarkonium helicity 
amplitudes are still rare on the market. In this presentation, we have achieved the first 
step to the development of an automatic Monte Carlo generator for heavy quarkonium. 
Our program is an extension of the present published HELAC [5j [6J, [T] , which is based on an 
off-shell recursive algorithm or Dyson-Schwinger equations. We dubed it as HELAC-Onia. 
It provides an automatic computation tool for heavy quarkonium helicity amplitudes in 
the standard model with high efficiency. We have also shown the applications of our tool 
to the various aspects of the heavy quarkonium production from pp,pp and e + e~ collisions. 

The following steps are to realize the automation of the next-to-leading order quarko- 
nium helicity amplitudes computations. With such a code, one can perform the analysis 
of the heavy quarkonium production at a full next-leading order level, which is much more 
reliable and useful especially at the LHC. 
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